Rational Construction of Stochastic Numerical Methods for Molecular Sampling

In this article, we focus on the sampling of the configurational Gibbs-Boltzmann distribution, that is, the calculation of averages of functions of the position coordinates of a molecular $N$-body system modelled at constant temperature. We show how a formal series expansion of the invariant measure of a Langevin dynamics numerical method can be obtained in a straightforward way using the Baker-Campbell-Hausdorff lemma. We then compare Langevin dynamics integrators in terms of their invariant distributions and demonstrate a superconvergence property (4th order accuracy where only 2nd order would be expected) of one method in the high friction limit; this method, moreover, can be reduced to a simple modification of the Euler-Maruyama method for Brownian dynamics involving a non-Markovian (coloured noise) random process. In the Brownian dynamics case, 2nd order accuracy of the invariant density is achieved. All methods considered are efficient for molecular applications (requiring one force evaluation per timestep) and of a simple form. In fully resolved (long run) molecular dynamics simulations, for our favoured method, we observe up to two orders of magnitude improvement in configurational sampling accuracy for given stepsize with no evident reduction in the size of the largest usable timestep compared to common alternative methods

Generating Generalized Distributions from Dynamical Simulation

We present a general molecular-dynamics simulation scheme, based on the Nose' thermostat, for sampling according to arbitrary phase space distributions. We formulate numerical methods based on both Nose'-Hoover and Nose'-Poincare' thermostats for two specific classes of distributions; namely, those that are functions of the system Hamiltonian and those for which position and momentum are statistically independent. As an example, we propose a generalized variable temperature distribution that designed to accelerate sampling in molecular systems.Comment: 10 pages, 3 figure

A molecular-dynamics algorithm for mixed hard-core/continuous potentials

We present a new molecular-dynamics algorithm for integrating the equations of motion for a system of particles interacting with mixed continuous/impulsive forces. This method, which we call Impulsive Verlet, is constructed using operator splitting techniques similar to those that have been used successfully to generate a variety molecular-dynamics integrators. In numerical experiments, the Impulsive Verlet method is shown to be superior to previous methods with respect to stability and energy conservation in long simulations.Comment: 18 pages, 6 postscript figures, uses rotate.st

Ab initio mass tensor molecular dynamics

Mass tensor molecular dynamics was first introduced by Bennett [J. Comput. Phys. 19, 267 (1975)] for efficient sampling of phase space through the use of generalized atomic masses. Here, we show how to apply this method to ab initio molecular dynamics simulations with minimal computational overhead. Test calculations on liquid water show a threefold reduction in computational effort without making the fixed geometry approximation. We also present a simple recipe for estimating the optimal atomic masses using only the first derivatives of the potential energy.Comment: 19 pages, 5 figure

On the long-time integration of stochastic gradient systems

This article addresses the weak convergence of numerical methods for Brownian dynamics. Typical analyses of numerical methods for stochastic differential equations focus on properties such as the weak order which estimates the asymptotic (stepsize h → 0) convergence behavior of the error of finite time averages. Recently it has been demonstrated, by study of Fokker-Planck operators, that a non-Markovian numerical method [Leimkuhler and Matthews, 2013] generates approximations in the long time limit with higher accuracy order (2nd order) than would be expected from its weak convergence analysis (finite-time averages are 1st order accurate). In this article we describe the transition from the transient to the steady-state regime of this numerical method by estimating the time-dependency of the coefficients in an asymptotic expansion for the weak error, demonstrating that the convergence to 2nd order is exponentially rapid in time. Moreover, we provide numerical tests of the theory, including comparisons of the efficiencies of the Euler-Maruyama method, the popular 2nd order Heun method, and the non-Markovian method

Active swarms on a sphere

Here we show that coupling to curvature has profound effects on collective motion in active systems, leading to patterns not observed in flat space. Biological examples of such active motion in curved environments are numerous: curvature and tissue folding are crucial during gastrulation, epithelial and endothelial cells move on constantly growing, curved crypts and vili in the gut, and the mammalian corneal epithelium grows in a steady-state vortex pattern. On the physics side, droplets coated with actively driven microtubule bundles show active nematic patterns. We study a model of self-propelled particles with polar alignment on a sphere. Hallmarks of these motion patterns are a polar vortex and a circulating band arising due to the incompatibility between spherical topology and uniform motion - a consequence of the hairy ball theorem. We present analytical results showing that frustration due to curvature leads to stable elastic distortions storing energy in the band.Comment: 5 pages, 4 figures plus Supporting Informatio

A Robust Numerical Method for Integration of Point-Vortex Trajectories in Two Dimensions

The venerable 2D point-vortex model plays an important role as a simplified version of many disparate physical systems, including superfluids, Bose-Einstein condensates, certain plasma configurations, and inviscid turbulence. This system is also a veritable mathematical playground, touching upon many different disciplines from topology to dynamic systems theory. Point-vortex dynamics are described by a relatively simple system of nonlinear ODEs which can easily be integrated numerically using an appropriate adaptive time stepping method. As the separation between a pair of vortices relative to all other inter-vortex length scales decreases, however, the computational time required diverges. Accuracy is usually the most discouraging casualty when trying to account for such vortex motion, though the varying energy of this ostensibly Hamiltonian system is a potentially more serious problem. We solve these problems by a series of coordinate transformations: We first transform to action-angle coordinates, which, to lowest order, treat the close pair as a single vortex amongst all others with an internal degree of freedom. We next, and most importantly, apply Lie transform perturbation theory to remove the higher-order correction terms in succession. The overall transformation drastically increases the numerical efficiency and ensures that the total energy remains constant to high accuracy.Comment: 21 pages, 4 figure